{\rtf1\ansi\ansicpg1252\deff0\deflang1033{\fonttbl{\f0\fswiss\fcharset0 Arial;}}
{\*\generator Msftedit 5.41.15.1507;}\viewkind4\uc1\pard\f0\fs20 % HM3nn0_Dmc.m neural network for three interferometer data (including the high mass end of the low mass data) by Kari Hodge\par
% last updated on June 21, 2007\par
% this time, replace all unknown parameters with 0\par
tic\par
\par
clear all\par
close all\par
\par
% load hM0 % the arrays that you are inputting are explained in highMinputformat.m and effective SNRsq.m\par
% load Eformatoutput Ps Pb\par
load Eformatoutput_diff_mc\par
\par
% find sizes of uploaded arrays\par
[NsP, Ns]=size(Ps); % [number of signal parameters, number of signal events]\par
[NbP, Nb]=size(Pb);\par
\par
if (NsP~=NbP)\par
    disp('error in input -- dimensions of signal parameters and background parameters must match')\par
end\par
\par
% concatenate the signal and background input for the nn\par
Input=[Ps Pb];\par
\par
% set the target vectors (signal's target is 1 and background's is 0)\par
Ts=ones(1,Ns);\par
Tb=zeros(1,Nb);\par
Target = [Ts Tb];\par
\par
%Divide the input and target samples into training, validation, and test samples\par
Q = Ns+Nb;\par
iitst = 2:4:Q;\par
iival = 4:4:Q;\par
iitrn = [1:4:Q 3:4:Q];\par
vali_P = Input(:, iival); vali_T = Target(:, iival);\par
test_P = Input(:, iitst); test_T = Target(:, iitst);\par
% eSNRfortest_ = eSNRsqSUM_(:, iitst)\par
train_P = Input(:, iitrn); train_T = Target(:, iitrn);\par
% only take out the parameters you want to look at:\par
vali.P = vali_P(2:NsP-1,:); vali.T = vali_T(1,:);\par
test.P= test_P(2:NsP-1,:); test.T = test_T(1,:);\par
ptrain = train_P(2:NsP-1,:); ttrain = train_T(1,:);\par
\par
% 'normalize' & preprocess the data to make it easier for the nn to understand\par
mm = [min(Input(2:NsP-1,:),[],2)  max(Input(2:NsP-1,:),[],2)];\par
% if (0)\par
%     [ptrain,meanp,stdp] = prestd(ptrain);\par
%     [ptrain,transMat] = prepca(ptrain,0.);\par
%     vali.P = trastd(vali.P,meanp,stdp);\par
%     vali.P = trapca(vali.P,transMat);\par
%     test.P = trastd(test.P,meanp,stdp);\par
%     test.P = trapca(test.P,transMat);\par
% end\par
\par
% find the correct dimensionality (first layer has same number of neurons as there are variables in the input:\par
Nv=length(Input(:,1))-2;\par
disp('hidden layer Nv+2 neurons')\par
% initialize the weights the same way each time\par
rand('twister',5489);\par
% create a feed-forward back propogation NN with log-sig layers:\par
net = newff(mm, [Nv,Nv+2,1], \{ 'logsig','logsig','logsig'\}, 'trainrp');\par
\par
%  train, validate, and test the nn\par
[net,tr] = train(net,ptrain,ttrain,[],[],vali,test);\par
\par
% plot the squared error vs epoch\par
figure(100)\par
plot(tr.epoch,tr.perf,tr.epoch,tr.vperf,tr.epoch,tr.tperf)\par
legend('Training','Validation','Test',1,'FontSize',18);\par
ylabel('Squared Error','FontSize',18); xlabel('Epoch','FontSize',18)\par
set(gca,'FontSize',14)\par
% now run the test data through the trained net:\par
tt = sim(net,test.P);\par
\par
% get the original signal and background from what the net found\par
isg = find(test.T > 0.5);\par
sg = tt(isg);\par
sgID = test_P(1,isg);\par
sgSNRe2 = test_P(NsP,isg);\par
ibk = find(test.T < 0.5);\par
bk = tt(ibk);\par
bkID = test_P(1,ibk);\par
bkSNRe2 = test_P(NbP,ibk);\par
\par
% find the few loudest background events\par
sortedbk = sort(bk, 'descend');\par
\par
% histogram NN output for signal and background\par
xlim = 0:0.05:1;\par
nsg = histc(sg,xlim);\par
nbk = histc(bk,xlim);\par
figure(101)\par
bar(xlim',nsg','k')\par
title('histogram of signal"s significance','FontSize',18)\par
set(gca,'FontSize',14)\par
figure(102)\par
bar(xlim',nbk','m')\par
title('histogram of background"s significance','FontSize',18)\par
set(gca,'FontSize',14)\par
\par
% apply cut on NN output, print out successes and failures\par
cut = 0.01;\par
css = length(find(sg > cut));\par
csf = length(find(sg <= cut));\par
cbs = length(find(bk <= cut));\par
cbf = length(find(bk > cut));\par
% [cbs cbf csf css]\par
disp([num2str(css) ' signal events with a significance greater than ' num2str(cut)])\par
disp([num2str(cbs) ' background events with a significance less than ' num2str(cut)])\par
disp([num2str(csf) ' signal events with a significance less than ' num2str(cut)])\par
disp([num2str(cbf) ' background events with a significance greater than ' num2str(cut)])\par
\par
%css/cbf\par
\par
% try the Discriminant analysis:\par
% class = classify(test.P',ptrain',ttrain','quadratic');\par
% Nc = hist(2*test.T+class',4);\par
% Nc\par
\par
% find which points were classified by the neural network(given a certain\par
% cut, loop over cuts)\par
% nns=find(tt>cut);\par
% nnstt=tt(find(tt>cut));\par
% nnb=find(tt<cut);\par
% compare with the original points\par
% oPs=find(test.T==1);\par
% oPb=find(test.T==0);\par
\par
toc\par
ROC1\par
\par
% extract all the parameters for the test points that the neural network\par
% classified, and the 21th row is the significance statistic:\par
sgArray=[test_P(:,isg);sg];\par
bkArray=[test_P(:,ibk);bk];\par
\par
[sortedsgArray,isort]=sort(sgArray(21,:));\par
sortedsgArray(1:24,1:2712)=sgArray(1:21,isort);\par
\par
[sortedbkArray,isortb]=sort(bkArray(21,:));\par
sortedbkArray(1:24,:)=bkArray(1:21,isortb);\par
}
 